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TECHNICAL MEMORANDUM 


THREE-DIMENSIONAL NUMERICAL MODELING OF MAGNETOHYDRODYNAMIC 

AUGMENTED PROPULSION EXPERIMENT 

1. INTRODUCTION 


Magnetohydrodynamic (MHD) augmentation of thermal propulsion systems has been suggested 
as a plausible means of boosting exhaust velocity, and possibly improving overall specific energy attri- 
butes. In this way, one might hope to reduce fuel fraction and shrink vehicle size without sacrificing 
payload delivery capability. To obtain a meaningful improvement in fuel fraction, however, it can be 
shown that the electrical augmentation power must be greater in magnitude than the thermal power of 
the unaugmented source. Such considerations lead to some extremely daunting technical challenges, 
particularly in relation to the development of an onboard electrical power source with adequate power 
density characteristics. Nevertheless, several technological avenues can be identified that someday may 
lead to innovative, compact, high-power electrical energy sources possessing the required attributes, and 
exploratory pursuit of fundamental technical feasibility is not without credible justification. Moreover, 
MHD accelerator technology has potential dual use application in ground-based hypersonic wind tunnel 
facilities where the electrical power source weight is of little or no concern. 

The essential requirement for using electromagnetic acceleration techniques is that the exhaust 
jet from the thermal propulsion source be electrically conductive. In practice, this can be accomplished 
by seeding the combustor flow of a chemical rocket with an alkali metal vapor, such as cesium, rubid- 
ium, potassium, and associated compounds. Because alkali metals have a relatively low ionization 
potential, the energy consumed in fully ionizing the seed is only a small fraction of the available thermal 
energy. Furthermore, the relatively low plasma working temperature is compatible with existing materi- 
als and regenerative cooling techniques. Using energetic rocket fuels, this method is known to produce 
supersonic plasma flows with an electrical conductivity on the order of 102 S/m, which is sufficient for 
evoking significant MHD interaction. At this level of MHD interaction, steady plasma acceleration is 
best invoked through externally imposed crossed electric and magnetic fields. This configuration gives 
rise to the so-called “crossed field MHD accelerator” in which the imposed Lorentz body force acceler- 
ates the flow. 

Small prototypes for this class of plasma accelerator have been designed and built, but almost 
exclusively from the standpoint of producing a hypersonic wind tunnel rather than a propulsive device. 
Testing with these prototype devices has clearly demonstrated flow acceleration, but diagnostic limita- 
tions have prevented complete delineation of the fundamental physical phenomena. Many uncertainties 
remain including the relative importance of electromagnetic versus electrothermal effects, achievable 
accelerator efficiencies, achievable current densities, maximum sustainable axial electric field without 
interelectrode arcing, effect of near-wall velocity overshoot phenomena, effect of microarcing in the cold 


electrode boundary layer, multiterminal loading of a segmented Faraday channel versus two- 
terminal loading of a diagonal wall configuration, and thermal loading and erosive effects with respect 
to long-term channel survivability, to name a few. 

Several years ago, NASA Marshall Space Flight Center (MSFC) initiated development of the 
Magnetohydrodynamic Augmented Propulsion Experiment (MAPX) for the purpose of resolving some 
of the critical technical issues associated with the use of MHD accelerators as thrust augmentation 
devices. A summary description of the project was previously published that included a thorough histori- 
cal account of preceding MHD accelerator research programs, a detailed description of the experiment 
configuration, and results from a preliminary quasi-one-dimensional engineering performance analysis 
and design study . 1 A recent status update for the project has summarized hardware design and develop- 
ment progress over the intervening time period . 2 In parallel with this work, a three-dimensional numeri- 
cal model based on the parabolized Navier-Stokes (PNS) approximation was developed to provide 
a comprehensive analysis of MHD accelerator performance. This new model was evolved from a preex- 
isting three-dimensional numerical model previously validated for MHD generator performance , 3 
and was modified to incorporate the NASA Chemical Equilibrium with Applications (CEA) code 4 as 
a means of calculating thermodynamic and species concentrations properties for the estimation of ther- 
moelectric properties of partially ionized gases using a fundamental technique based on electron-neutral 
momentum transfer cross sections . 5 The purpose of this Technical Memorandum (TM) is to summarize 
development of the parabolic three-dimensional numerical model, present pretest analyses and perfor- 
mance optimization results, and provide a recommended test configuration for initial MAPX tests. 
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2. EXPERIMENT DESCRIPTION 


The MAPX facility adapts a traditional linear MHD flow path configuration, as shown schemati- 
cally in figure 1. First, the working fluid (i.e., nitrogen) is heated in a 1.5-MW e segmented multigas arc- 
heater to a stagnation temperature, Tq~ 4,000-4,500 K at a stagnation pressure, P () ~ 1 0 atm. The hot 
gas then enters a mixing chamber where alkali metal seed (i.e., NaK) is injected into the flow stream 
after which it is expanded through a primary nozzle to a Mach number in the range of M~ 1 .25 to 1.5. 

A 2-MW e MHD accelerator directly increases the energy and momentum of the flow, which is further 
diffused in a secondary nozzle to obtain the maximum possible jet velocity. The secondary nozzle 
exhausts into a large windowed test section equipped with a stinger mounted stagnation probe or aero 
model. The test section is attached to a nitrogen-driven ejector pump designed to maintain a backpres- 
sure in the range of 1 to 3 psia. 


2-MW e MHD Accelerator 
and 2-T Electromagnet 



Figure 1 . Schematic of the NASA MSFC MAPX facility. The major flow path components are 
(1) 1.5-MW e arc-heater, (2) seed injector and mixing chamber, (3) primary expansion 
nozzle, (4) 2-MW e MHD accelerator channel and 2-T electromagnet, (5) secondary 
nozzle, (6) windowed test section, and (7) nitrogen-driven ejector pump. 


2.1 Hot Gas Source and Entrance Path 

The hyperthermal stagnation conditions are generated by a 1.5-MW e (nominal) segmented multi- 
gas arc -heater, which operates in a wall-stabilized, constricted arc, dc discharge mode. The arc-heater 
is energized by a saturable reactor dc power supply that can sustain a continuous operating power of 
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0.75 MW e on an indefinite basis and can deliver an intermittent power burst of 1 .5 MW e for 5 to 10 min. 
The purpose of the thermal driver is to maximize the MHD interaction by delivering the highest possible 
mass throughput and maximizing the accelerator channel size. 

This particular arc-heater has a 1-in internal bore diameter and follows the traditional segmenta- 
tion design philosophy, whereby alternating conductor/insulator wafers are stacked together to form the 
full length assembly, as illustrated in figure 2. The 3/8-in-thick, heat-conducting copper segments are 
water cooled and are separated by boron nitride insulators in stacked pack subassemblies, which are held 
securely together by four INCONEL® tie rods. These subassembly packs are then attached in a sequen- 
tial manner to form the full arc-heater column, which spans an overall length of =1 m in the three-pack 
configuration shown. The working gas is injected tangentially through four 0.048-in jets in a primary gas 
injection segment near the rear of the arc-heater, and a dc arc discharge is established between a tungsten 
cathode button in the rear sealing flange and a copper anode ring at the arc-heater exhaust. A magnetic 
spin coil is located around the anode ring to induce continuous rotation of the arc attachment point. 


Primary Gas 
Injection 
Tungsten Segment 
Cathode 
Button 


Starting Segment 


Vacuum 
Contactor 

— OO — V\AA 


| Secondary 
Gas Injection 
Segment 


3 0) 


Spin 

p^ oil oAnode 
1 — J-’Ring 



-V 2,500-5,000 V 


Saturable Reactor Power Supply 
0.75 MW e (Continuous) 

1.5 MW e (5-10 min) 



Typical Segmentation Assembly 


Figure 2. Segmentation and assembly detail of the MSFC l-MW e multigas arc -heater. 


An exploded schematic of the entrance flow path assembly is shown in figure 3. This assembly 
accomplishes the following primary functions: (1) Electrical isolation of the arc -heater from the MHD 
accelerator and circular-to-square flow path transition, (2) NaK seed injection and mixing, and (3) flow 
acceleration via the nozzle. The isolation flange is designed for heat sink operation and is an expendable 
item, with the other components being water-cooled copper pieces intended for extended service. An 
alumina insert is located inside the isolation flange (as shown in fig. 3) that provides the geometric tran- 
sition of the flow passage from the circular geometry of the arc-heater exit to the rectangular geometry 
of the MHD accelerator. This geometric transition is facilitated by the internal loft (flow area distribu- 
tion) of the alumina insert that gradually transforms from circular to square. It is also worth noting that, 
because the circular-to-square transition occurs along the length of the isolation flange assembly without 
a change in diameter (width), the cross-sectional area increases by a factor of 4 /%. 
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Primary Nozzle 
Assembly 


Arc-h 



Figure 3. Exploded schematic of MAPX entrance flow path assembly. 


Over the years, a considerable amount of historical performance data have been accumulated 
and cataloged for this particular arc-heater, and the resulting database can be used to project perfor- 
mance characteristics at representative MAPX test conditions. Projected variation in electric-to-thermal 
conversion efficiency with applied electrical power, for instance, is shown in figure 4, along with some 
recently acquired experimental data. These experimental efficiencies were inferred from calorimeter 
measurements obtained during high flow rate nitrogen runs using an uncooled graphite nozzle with 
a 5/8-in 0 throat. In general, the projected efficiencies were within one to two percentage points of 
the experimentally observed values over the examined power range. 

Preliminary analysis, assuming an applied electrical power of 1 . 1 MW e and nozzle losses of no 
more than 20% of the available total enthalpy, indicated that satisfactory conditions could be achieved 
with 130 g/s of nitrogen and 1.5% NaK (by weight) using a 0.567x0.567 in 2 throat with an area expan- 
sion ratio of A/A* = 1.142297. The 20% nozzle heat loss limit was derived from practical experience 
with water-cooled copper nozzles in high-temperature combustors. The primary nozzle performance 
characteristics were estimated using a modified version of the CEA code, which incorporates a method 
for computing plasma electrical transport properties. At the 1.1 -MW e operating design point, the arc- 
heater electric-to-thermal conversion efficiency was experimentally determined to be =63%, yielding 
an equilibrium total temperature of 3,302 K at the primary nozzle entrance. Anticipated thermodynamic/ 
electrical conditions at the accelerator entrance are summarized in figure 5 as a function of actual pri- 
mary nozzle performance. These chemical equilibrium calculations indicate that acceptable accelerator 
inlet conditions can be attained even for the worst-case performance scenario. In this performance- 
limiting situation, the inlet static temperature and electrical conductivity are 2,397 K and 1 1 S/m, 
respectively. 
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Figure 4. Electric-to-thermal conversion efficiency with assumed nozzle cooling loses. 
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Figure 5. Primary nozzle performance characteristics. 
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2.2 Magnetohydrodynamic Accelerator Channel 
2.2.1 Electrical Loading Configuration 

Alternative configurations for linear MHD accelerator channels are depicted in figure 6 where the 
optimal MHD accelerator configuration is determined by the ultimate application needs. From a perfor- 
mance standpoint, the Hall configuration (fig. 6(a)) is more effective for low-density flows whereas the 
Faraday configuration (fig. 6(b)), with segmentation to neutralize the Hall current, is superior for high- 
density flows. The major drawback of the Faraday configuration, however, is the separate power condi- 
tioning required for each electrode pair which leads to a complex and expensive system. In many cases, 
particularly flight applications, multiterminal loading is not practical. 

Alternative two-terminal loading schemes have been proposed to avoid the multiterminal com- 
plications while attempting to reap the major benefit associated with the Faraday configuration; i.e., Hall 
current neutralization. For example, the standard segmented Faraday channel may be externally diago- 
nalized in a series connected scheme (fig. 6(c)), or one could adopt a diagonal conducting wall (DCW) 
configuration in which slanted window frame-like electrode elements are stacked with thin insulators 
to form a complete channel (fig. 6(d)). The DCW configuration not only simplifies fabrication and 
improves strength but provides superior performance to the externally shorted (i.e., series connected) 
device by allowing current to flow to the sidewalls. 



A A A A A A A 



vvvvv 
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(a) Linear Hall Accelerator 


(b) Segmented Faraday Accelerator 



(c) Series Connected Diagonal Accelerator 



uxB 

(d) Diagonal Conducting Wall Accelerator 


Figure 6. Alternative design configuration for linear MHD accelerator channels: (a) Linear 
hall accelerator, (b) segmented Faraday accelerator, (c) series connected diagonal 
accelerator, and (d) diagonal conducting wall accelerator. 
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Although the DCW configuration is the overall best candidate for flight implementation, 
the MAPX design was based on an externally diagonalized series connected configuration for reasons 
of cost and flexibility (e.g., effective wall angle adjustability). 

2.2.2 Engineering Design and Performance Analysis 

Numerous investigations have clearly established that MHD channel flows are subject to sig- 
nificant three-dimensional effects. Thus, averaging the governing MHD equations (magnetic Reynolds 
number «1) to obtain a quasi-one-dimensional engineering model requires the adoption of major 
simplifying assumptions. Nevertheless, many of these assumptions including wall friction, wall heat 
flux, and near-electrode voltage drops may be approximately accounted for through the introduction 
of appropriate physical wall functions for the boundary layer. 6 ’ 7 

Engineering design and performance analyses of the MAPX accelerator were carried out using 
a legacy code based on one such approach. This engineering code was initially developed within the 
Energy Conversion Division at the University of Tennessee Space Institute (UTSI) in support of the 
Department of Energy MHD Commercial Power program. Over the years, the code was evolved and 
expanded to encompass a range of generator and accelerator loading configurations. In most respects, 
the development is similar to that described for other nonperfect-gas, quasi-one-dimensional analyses, 
the principle idiosyncrasies being associated with the physical submodeling. The code solves the gov- 
erning internal duct flow equations for conservation of mass, momentum, and energy together with the 
equation of state and boundary layer wall functions using a fourth-order Runge-Kutta numerical inte- 
gration scheme. It uses a real-gas equation of state and assumes local thermodynamic equilibrium as 
predicted by the CEA code with appropriate modifications for computing electrical transport properties. 
Values for the empirical constants associated with various physical submodels have been established 
through extensive benchmarking experience. 

Heat transfer and frictional wall losses are computed intrinsically in the code and require input 
of the wall temperature and roughness height. Near-wall electrical losses are also treated intrinsically 
through integration of the conductivity profile as defined by velocity and thermal boundary layer cor- 
relations for fully turbulent flow. This correlation computes the boundary layer growth along the MHD 
accelerator duct through definition/input of the initial boundary layer height and shear (viscosity as 
a function of temperature). The velocity and temperature profiles are taken as 1 In power-law distribu- 
tions. This approach also relies on user specification of the Rosa G factor to account for plasma nonuni- 
formities and effective voltage drop. A value of G ~ 2 was anticipated based on past experience and was 
therefore utilized for engineering design and performance analyses purposes. 

Actual accelerator performance depends, of course, on constraints imposed by the available 
magnet and power supply equipment. To meet the research goals of this program, a water-cooled, 2-T 
electromagnet was acquired from UTSI and refurbished to support general MHD research at MSFC. 

The specifications for this magnet are summarized in table 1 . A new 3, 000- A, 75-V dc power supply 
was acquired to power the magnet, and the entire system has been installed and integrated into the 
MAPX flow train. A shakedown test of the electromagnet system was conducted and the measured field 
profile at 2,400 A is shown in figure 7. A 2-MW e , high-voltage dc power supply was also acquired and 
installed to power the accelerator. The voltage on the unit is variable to 10 kV and is capable of deliver- 
ing 300 A at 6,700 V. 
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Table 1. Electromagnetic performance specifications. 


Magnetic field strength 

2 T 

Air gap 

4 in 

Pole cap length 

36 in 

Voltage 

65 V 

Maximum current 

2,400 A 

Cooling water (70 psig supply) 

50 gal/min 



Axial Position (in) 

Figure 7. Measured centerline magnetic field strength profiles at 2,400-A applied current. 

The detailed design process entailed several iterative calculations in an attempt to optimize 
stagnation pressure rise by varying load current, channel divergence, and channel length. These engi- 
neering calculations were carried out utilizing the measured centerline magnetic field profile at 2,400 A 
The height-to- width aspect ratio was unity at the inlet, and the E-field and B -field walls were diverged 
to accommodate boundary layer growth and flow expansion. Preliminary analysis of the axial current 
neutralized operation indicated that the optimal diagonalization angle was near 6 =-45°. The resulting 
physical specifications selected for the accelerator design based on this extensive engineering study are 
summarized in table 2. 
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Table 2. MAPX accelerator specifications. 


Inlet heightxwidth 

1.6x1. 6 cm 2 

Channel divergence 

1.0° 

Electrode width 

1.0 cm 

Insulator width 

0.5 cm 

Active length 

90 cm 

Powered electrodes 

60 

Total length 

96 cm 

Total electrodes 

65 

Exit heightxwidth 

3.6x3.6 cm 2 

Seed (NaK) 

1.5% 

Nitrogen flow rate 

130 g/s 

^0,in 

3,120 K 

^0,in 

8.5 atm 

"in 

1,312 m/s 

^in 

25 S/m 

An 

0.7 
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3. MODEL DESCRIPTION 


It is widely recognized that the current and electric field structure in MHD devices is inherently 
three-dimensional and spatial separation of physical processes is not warranted in general. Generally, 
it is sufficient to invoke an infinite segmentation assumption (i.e., the streamwise variation in electri- 
cal properties is small in comparison to transverse variations) and couple an approximate cross-plane 
electrical model with a three-dimensional flow analysis, as exemplified for generator configurations by 
Bityurin et al. 8 and Ahluwalia et al. 9 

Here, we adopt the multigrid magnetohydrodynamic (MGMHD) model 3 previously developed 
by Argonne National Laboratory for the analysis of MHD generators and adapt it for the MHD accelera- 
tor configurations. The resulting development, to be described in this TM, is designated as the universal 
MGMHD model (UMM). 10 UMM is a PNS computational fluid dynamics- (CFD-) based methodology 
capable of simulating both generator and accelerator modes of operation in Hall, Faraday, and diagonal 
configurations. The numerical model utilizes the CEA code for calculating thermodynamic properties 
and species concentrations and adopts a fundamental kinetics model for the estimation of electrical 
transport properties from experimental electron-neutral momentum transfer cross sections. Section 3 
provides a brief overview of the MGMHD numerical model, including a detailed discussion of its capa- 
bilities and limitations. 


3.1 Baseline Modeling Capabilities 

The MGMHD computer code is an updated version of the single-grid, three-dimensional magne- 
tohydrodynamic (TDMHD) code developed at the Argonne National Laboratory in 1982 for the analysis 
of MHD generators and diffusers. 3 The MGMHD code retains the TDMHD formulation of three- 
dimensional partial differential equations for flow and electrical fields but incorporates an advanced 
multigrid solution algorithm to reduce computer convergence time. The MHD equation set comprises 
the mass continuity equation, three momentum equations, the energy equation, two turbulence model 
equations, and Maxwell’s equations. Turbulence is represented by a two-equation model ( k - s ) in which 
partial differential equations are solved for the turbulence energy (k) and its dissipation rate (e). Lastly, 
the MGMHD code provides two unique features: (1) A full approximation storage (FAS), block implicit 
multigrid, finite-difference solution procedure for the cross-stream hydrodynamic equations and 
(2) a FAS multigrid, finite-difference solution procedure for cross-stream electrical potential equations. 

3.1.1 Fluid Dynamic Model 

Flow processes in an MHD generator can be generally represented by the three-dimensional 
Navier-Stokes equations. In an MHD channel, however, the flow is predominantly in the stream-wise 
direction, which allows for certain simplifications when considering the order of magnitude of various 
terms. The simplification made here is referred to as the PNS approximation, which is akin to the 
boundary layer approximation made in solving problems such as the flat-plate flow, and consists of the 
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following major assumptions: (1) Diffusion flux in the stream-wise coordinate may be neglected and 
(2) the pressure gradient in the axial momentum equation is assumed to be uniform over the channel 
cross section. 

The PNS approximation permits the solution to be marched forward plane-by-plane from the 
inlet to the exit of the channel. This approximation does not introduce any significant errors when the 
flow is predominantly in one direction, and it eliminates the need to iterate between the inlet and exit of 
the channel, thereby making the calculation converge more quickly. In this case, the pressure field, P, is 
split into two components: (1) Pi oca \, the local cross-sectional pressure, and (2) p avg , the cross-sectional 
average pressure, which, as the following equation shows, is only a function of the axial coordinate: 

P(x,y,z) = P local (x,y,z) + p avg (x) . (1) 

The coordinate system utilized by the MGMHD numerical model is Cartesian; wherein, the mass 
continuity, momentum, and energy (enthalpy) equations are represented by the following equation set: 

Mass: 


sH + £h + !m=° 


x-Momentum (axial direction in channel): 


j’-Momentum: 


z-Momentum: 


3/ \ d j \ d i \ 3/> 3 f 3 t 


dx 


d_ 

dx 


(puv^+-^-(pwJ+-^-(pwvj =— ^i° cal -l yy -l yz 


d_ 
dz ' 


dr dT, 


dy dy dz y 


3/ \ 3 / \ 3 / \ ^Pineal zy 3 T 

— (puw) + — ( pvw ) + — [pww) = locai_ + y_ + — ZZ 

1 ’ dy { ’ dz { ’ dz dy dz 


dx 


Enthalpy (energy): 


( 2 ) 


(3) 


(4) 


(5) 


± ipuh) + ± (pvh) + l (pwh) = -^-^ + u ^ + v ^ + w ^L + ^ + D 

dx dy dz dy dz dx dy dz O 

where the viscous dissipation rate, D, is calculated as follows: 


( 6 ) 


D = (v l + H t ) 
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In the preceding equations, u, v, and w are the x, y, and z components of the u velocity vector, 
respectively. The symbol p represents mass density, r represents shear stress, a is electrical conductivity, 
h is enthalpy, J is current density, and B is the magnetic field (here, B represents only the scalar mag- 
nitude of the magnetic field). The MGMHD code assumes steady state; therefore, any time-dependent 
fluctuations will not be captured with this numerical model. Furthermore, note that equation (6) uses 
enthalpy instead of c p T , the product of specific heat and temperature, in the energy equation. This is 
because the MGMHD code advances the energy equation using enthalpy, not temperature. Advancing 
the energy equation using temperature requires that one assume a constant specific heat (i.e., taking c p T 
out of the derivative in energy equation), which is not true in some temperature ranges. Advancing the 
energy equation using enthalpy (or internal energy) requires no such assumption. 

3.1.2 Electrical Model 

Slow calculation of the electrical potential solution is one difficulty in performing a complete 
three-dimensional MHD analyses with finite channel segmentation. In this section, the cross-sectional 
infinite segmentation electrical model, developed by Ahluwalia 9 and utilized in the MGMHD model, is 
summarized. In the infinite segmentation model, axial variations of the flow and electrical variables are 
negligible in comparison to their cross-plane variations. Thus, a solution may be obtained by marching 
stepwise from the entrance of the channel to the exit. This model was successfully applied to the three- 
dimensional analysis of Faraday, diagonal insulating, and diagonal conducting sidewall MHD generators 
and carefully validated against experimental data with good success. The electrical model is coupled to 
the fluid dynamic model presented above. More extensive discussion and detailed derivations may be 
found elsewhere. 9 ’ 10 

In the MGMHD code, the electrical governing equations consist of Maxwell’s steady-state equa- 
tions and Ohm’s Law. The electrical field, E, and the current density, J, are obtained by solving the fol- 
lowing equations: 


VxE = 0 , 

(8) 

V-J=0 , 
and 

(9) 

J = ct(E-i-uxB)-— JxB . 
v ’ B 

(10) 

In the above equation, /? is the Hall parameter. As typical for MHD channel flows, these equa- 
tions are valid for low magnetic Reynolds number flows, where the induced magnetic can be neglected. 
Equations (8) and (9) make it possible to define an electric potential, y/, such as 

> 

i 

ii 

W 

(11) 


which, because of the steady-state assumption, effectively reduces the vector problem to that of a scalar. 
Also, note that the MGMHD code assumes that the magnetic field, B, is sectionally uniform and oriented 
in the positive z-direction, as shown in figures 8 and 9. 
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Electrode 


Electrode 



Figure 8. Schematic representation of MHD generator configurations in Cartesian coordinates: 
(a) Diagonal conducting wall, (b) diagonal insulating sidewall, and (c) Faraday. 


The infinite segmentation model is constructed from equations (8) and (9) by neglecting axial 
variations of the flow and electrical variables in comparison with their cross-plane variations. Conse- 
quently, the axial electric field, E x , is assumed to be constant in a cross-sectional plane. By defining 
a set of oblique coordinates in the general case of a diagonal connection, the electric potential can 
therefore be represented by the relation 

yf = -(x'-ycot6 d )E x + <p(y,z) , (12) 

where 6 d is the diagonalization angle (measured in a counterclockwise direction from the positive 
x-axis), and cp(y,z) is a two-dimensional potential. The first term on the right side of equation (12) rep- 
resents the contribution of the axial (x-axis) electric field. Substituting equation (12) into equations (8)- 
(10) results in the following equation for the two-dimensional function, <p(y,z ): 


d_ 

dy 


d(p 

r " 0T 


+ • 


dz 


d(p 

dz 


=£ *|;[(^" cotS d CT J-|;[ <T » s (*'-£ v )] • 
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where (7 w =(t/(1+/3 2 ). Since E x is not known a’ priori, and because equation (13) is linear in (p, 
a decomposition of the potential is permitted as follows: 


(p(y,z) = E x (p l +(p 2 


(14) 


The two functions, (p { and (p 2 , are governed by the following equation set: 

i W=|;[(0- co<( O <7 »] 

and 


where the operator, L, is 
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f d) 
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dy 
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dz 

i dz J 


(15) 


(16) 


(17) 


Note that the function cp l is driven by temperature nonuniformities, while <p 2 is driven by temper- 
ature and velocity nonuniformities. Both functions are independent of E x , which is reduced to a param- 
eter of the problem to be determined later from the specification of the external electrical connection and 
loading. From Ohm’s Law, the current density components — J x , J y , and J z — can be expressed as 


15 



( 18 ) 



(19) 


and 



( 20 ) 


To find the cross-sectional averages of J x and J y , the cross-sectional spatial average must be defined as 


where A is the cross-sectional area normal to the channel, which yields the following average equations: 


Applying the boundary conditions for diagonal insulating sidewall and diagonal conducting sidewall 
MHD generators, the axial electric field is 


3.1.3 State Equations 

In order to close the system of fluid dynamic and electrical equations, it is necessary to establish 
a relationship between the fluid properties and the chosen solution transport variables ( u , v, w,p avg , h ). 
In the MGMHD code, electrical and thermodynamic properties are calculated from two independent 
thermodynamic variables, (p avg , h) or {p avg , T), depending on which of the four following options are 
chosen by the user: 



( 21 ) 



( 22 ) 


and 



(23) 


E 



( cr w (/ j + cot ^)( < P 2 J + ^))~- 5 ( (7 H( 1 -/ 3cot ^) v ) 

( 1 + cot 2 3 0^ ) (<7 M (/3 + cot 9 d ) <p, v ) 


(24) 
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(1) Constant fluid properties; 

(2) Variable fluid properties with constant electrical properties; 

(3) Variable fluid properties with o and /3 calculated from empirical correlations; and, 


(4) Variable fluid and electrical properties calculated from fourth-order interpolation 
polynomials in p avg and h. 

Each of the above-stated property options is detailed in the MGMHD user manual; 3 however, it is 
important to note that all of the options use approximations and/or polynomial curve fits to generate the 
thermodynamic and electrical properties. The baseline MGMHD code does not include capability for 
generalized calculation of electrical transport properties in partially ionized gases. 

3.1.4 Required Modifications and Improvements 

The baseline MGMHD code is inadequate for three-dimensional modeling of MHD accelerators; 
however, it does fulfill some of the basic requirements and its structure and flexibility allow for the 
necessary modifications. As noted above, the PNS formulation of the MGMHD model enables three- 
dimensional calculations to proceed without the penalty of slow convergence times or the need for 
massive computing resources. Moreover, the MGMHD code is mathematically structured for diago- 
nal configurations, although the diagonalization angle must remain constant over the entire length of 
the MHD channel. Unfortunately, the thermodynamic and electrical conductivity submodels in the 
MGMHD code are, at best, generic curve fits, based on a few user-defined constants, but the code does 
allow the user to supply their own thermodynamic and electrical transport data via an input file. 

In order to properly model MHD accelerators, the MGMHD code requires certain specific 
modifications, as follows: 

• Structural modifications to permit both generator and accelerator calculations in Faraday, Hall, 
and diagonal configurations. 

• Modifications to allow arbitrary variations in diagonalization angle along the channel length. 

• Addition of thermodynamic and electrical transport property subroutines based on a generalized 
kinetics model applicable to arbitrary partially-ionized gas mixtures. 

Though not necessarily required for the modeling of diagonal accelerators, there are several other 
desirable features that would enhance and improve the capabilities of the MGMHD code. These addi- 
tional improvements are as follows: 

• The addition of a realistic power-takeoff scheme at the entrance and exit of the simulated MHD 
channel which would allow for tailored distribution of current flow in these regions help and thereby 
avoid MHD compression effects at the inlet. 

• Incorporation of accelerator efficiency calculations which would allow the user to more quickly 
determine optimum accelerator configurations. 

• A three-dimensional calculation of total pressure as an aid to design optimization. 
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An axial current neutralized optimization option (an extension of the varying diagonal angle), which 
would facilitate determination of an optimum diagonalization angle and minimization of axial ( J x ) 
current flow. 


3.2 Extended Modeling Capabilities 

As noted above, the baseline MGMHD code was inadequate for proper modeling of diagonalized 
MHD accelerators, and significant modifications were required to correct these deficiencies. The struc- 
ture of the modified numerical model, applicable to both MHD generators and accelerators in Faraday, 
Hall, and diagonal configurations, is referred to as the universal MGMHD model (UMM) and may be 
separated into two major blocks, as shown in figure 10. The thermoelectric block consists of the thermo- 
dynamic submodel (i.e., CEA routine), the electrical transport property submodel, along with library 
and input files that control the overall thermoelectric computation process. The MHD block consists 
of a modified version of the baseline MGMHD code. 


UMM Coding Blocks 


MHD 


Thermoelectric 


mhncom.h 
angle. in 


( \ 

i 




MHD.f 


THERMHD.f 


CEB.f 


OMEGA.f 


ECON.h 

mhd.in 

thermhd.in 

mapx.in 




thermo. inp 
trans.inp 


Figure 10. Structural coding blocks for the UMM three-dimensional numerical model. 


3.2.1 Thermodynamic and Electrical Transport Property Models 

The program thermhd.f manages the three routines named ceb.f, omega.f, and econ.h as shown in 
figure 10, where the “.f” extension denotes a Fortran file and the “.h” denotes an include file. The input 
file for thermhd.f is thermhd.in, in which the user defines the pressure and enthalpy ranges of interest 
and the desired tabulation intervals. The primary purpose of the thermoelectric codes is to create a look- 
up table of thermodynamic and electrical transport property data (temperature, molecular weight, ratio 
of specific heats, laminar viscosity, electrical conductivity, and total charge) that can be accessed by the 
MHD numerical block as it solves the coupled fluid dynamic and electromagnetic equations. The ceb.f 
routine has three input files: (1) thermo. inp, which automatically creates a library of thermodynamic data 
(no user input required), (2) trans.inp, which automatically creates a library of transport property data 
(also, no user input necessary), and (3) mapx.in, where the user defines the fluid species, weight percent- 
ages, and initial temperatures. The major difference between the ceb.f routine and the CEA code is that 
ceb.f has been modified to output an unformatted plot file of 39 variables representing the 9 thermody- 
namic variables and 30 gas species concentrations for which electron-momentum, cross-sectional data 
are available. It is important to note that all changes made to the CEA code deal with input/ouput (I/O) 
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and formatting (i.e., none of the thermodynamic equations in the CEA code were altered). The ceb.f 
code runs at every step in the pressure and enthalpy looping process, as described in subsequent sec- 
tions. For every step in the pressure and enthalpy loops, the ceb.f code provides the following informa- 
tion to the MHD numerical model: 

• Temperature (7). 

• Molecular weight (MW). 

• Ratio of specific heats (y). 

• Laminar viscosity (/q). 

• Total charge (eN e ). 

In the UMM code, electrical conductivity is calculated in two parts, based on species colli- 
sions, using the omega.f and econ.h routines. Specifically, omega.f is a self-contained (i.e., no input file 
required) Fortran file that calculates electron-neutral collision integrals for 30 specified species using 
momentum-transfer collision cross-section data. 11 ’ 12 Charged-particle collision integrals are then 
estimated in econ.h. All the collision integrals along with the species and total number densities men- 
tioned earlier are then used by econ.h, which is an include file called from thermhd.f, to calculate the 
electrical conductivity of the gas mixture, based on a method described by Devoto, 5 from the fundamen- 
tal equation 


(7 = 


e 2 n 


pkT M 


'5'n.m.Z D . 

.i j j i. / 


(25) 


where e is the electron charge, n is the total number density, p is the total mass density, rij is the number 
density of the yth species, Dy is the multicomponent diffusion coefficient, and j is the counter, with 1 
representing electrons and 2 through C being ions (where >2 denoted multiple ionic states). 


The transport properties and diffusion coefficients in the UMM model are derived from the 
collisional Boltzmann equation: 


df F df 

dt m d\ 


\ dt Jc 


(26) 


One solution to the Boltzmann equation is the Chapman-Enskog perturbation method (which is rigor- 
ously explained and expanded in Hirschfelder, Curtis, and Bird 13 and Chapman and Cowling 14 ) where 
one assumes that collisions are the driving mechanism for translational nonequilibrium of the gases and 
are such that the distribution functions of the various species differ only slightly from a Maxwellian 
distribution function. The Chapman-Enskog method uses a perturbation factor based on a scale length of 
macroscopic quantities and the mean free path between collisions. The distribution function is expanded, 
order by order, in the small perturbation parameter: 14 


/(r,u,t) = / 0 (r,u,f) + e/j (r,u,t) + £ 2 / 2 (r,u,t) ■ 


(27) 


where / 0 ,/j, and / 2 are assumed to be of the same order of magnitude. If there are no gradients in the 
composition, velocity, and temperature in the gas, this reduces to a Maxwellian distribution: 
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/ 0 (r,u,t) = rc(r) 
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(28) 


Typically, one is interested in the properties of gases that are under conditions only slightly dif- 
ferent from equilibrium. In fact, it is only under these conditions that the flux vectors are linear in the 
derivatives and the usual definitions of the transport coefficients apply. 13 ’ 14 In this limit, the distribu- 
tion function is nearly Maxwellian, and the Boltzmann equation can be solved by the Chapman-Enskog 
perturbation method. The resulting solutions are then used to obtain expressions for the fluxes and the 
transport coefficients. 13 ’ 14 


For MHD channel flow, the following assumptions may be applied: 

• Inelastic collisions are negligible. 

• The number of ions and electrons are equal. 

• There is only one user-defined ion species, which is singly ionized. 

• The ion and electron temperatures are equal. 

• The electron-heavy (i.e., electron-neutral, electron-ion) collision terms are neglected in deriving the 
expression for the ion and atom transport properties (i.e., electrons do not alter the direction, velocity, 
or momentum of heavy particles). 

• The change in the heavy particle perturbation term during a collision will be neglected in obtaining 
the expressions for the electron transport properties (i.e., when an electron impacts a heavy particle, 
the perturbation to the heavy particle is insignificant when compared to the electron perturbation). 

• The contribution of ions to the electrical conductivity equation has been neglected (which is consistent 
with the simplification introduced in solving the linearized electron Boltzmann equation). 


Under these assumptions, the electron-neutral conductivity equation reduces to the form 
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(29) 


The solution of equation (27) for jqfr, v, t) is accomplished through an expansion of the func- 
tion in a finite series of Sonine polynomials, 5 with the level of approximation of the transport coeffi- 
cients being the number of terms in the expansion series. It has been shown that at least the second-order 
approximation is required for a reasonable level of accuracy. (See ref. 15.) In the fourth-order approxi- 
mation, the ordinary diffusion coefficient from equation (29) is given by 
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where \q\ is the determinant of the q mn elements. The q mn elements are the coefficients to the simpli- 
fied electron equilibrium Maxwellian distribution function and are dependent on Qy ’ s ^ , the Maxwellian 
velocity-averaged collision cross section between species i and j. For example, one of the more simple 
equations for q mn is given by 
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(31) 


where the equations for q mn become more complex as the values of m and n increase. The Maxwellian 
velocity-averaged collision cross section can be calculated through the relation 
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(32) 


Here, $ 2 =n ij g 2 /2kT, where Hy is the reduced mass of the colliding species i and j, and g is the initial rela- 
tive velocity before impact. Fortunately, it turns out that all of the Qy’^ {t ] terms from equation (32) 
with /> 1 have coefficients proportional to a power of m x hn - (where 1 = electrons and can be 
neglected. 16 This leaves only values of Ql., which are the well-known momentum-transfer collision 
cross sections, which can be experimentally measured as a function of incident electron energy . 11 ’ 12 
The integral in equation (32) takes the form of a half-range Gauss-Hermite polynomial, which is a 
special case of a Gauss-Laguerre integral equation, and can be solved using a 16-point Gauss-Laguerre 
quadrature (which is done in omega.f), essentially calculating the values of the average collision cross 
section for the selected species. From these cross-sectional values, the electrical conductivity based on 
electron-neutral collisions can be calculated from equations (29)-(31). It is known that in some gases, 
the momentum- transfer approximation is not accurate to the 1% level and the angular distributions must 
be taken into account. For MHD applications, however, the use of momentum-transfer cross sections is 
adequate. 


Because the ionization level in MHD devices is so low (say, on the order of l%-2%), electron- 
neutral collisions obviously play the dominant role in determining the electrical conductivity of the fluid. 
However, for completeness, the UMM code has a conductivity model that incorporates charged-particle 
collisions into the global electrical conductivity of the flow. Charged-particle collision integrals are cal- 
culated using a method described by Zollweg, 17 who modified the Spitzer formula to permit its applica- 
tion in calculating electrical conductivities of partially ionized gases in the ideal and nonideal plasma 
regions. Zollweg’s modifications to the Spitzer formula result in the relation 

Qy = 67r6 ln|l + 1.4A^ j 1/2 , (33) 

which approximates the momentum-transfer cross section for charged-particle collisions. This result 
is then used in conjunction with equations (32) and (29) to determine the electrical conductivity due to 
charged-particle collisions. Note that b is the averaged impact parameter and A is the Coulombic loga- 
rithmic term. 
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3.2.2 Electrical Conductivity Sequence 

Electrical conductivity is calculated in the UMM code using the omega.f and econ.h routines. 
Specifically, equation (25) is used to sum the different contributions from the different species interac- 
tions. Equation (32) calculates a collision integral using two methods, depending on the species involved 
in the collision. Electron-neutral collisions are calculated using a 16-point Gauss-Laguerre quadrature, 
where the momentum-transfer cross-section data comes from experimental results for the 30 available 
species. 11 ’ 12 Charged-particle collisions are approximated using the method devised by Zollweg, 17 
which is based on a modified Spitzer formula and assumes charge neutrality. 

The calculation of the integral cross section for the charged-particle interactions is fairly straight- 
forward where each cross section is calculated at the temperature of interest. However, the execution of 
the Gauss-Laguerre quadrature in omega.f is slightly more complex. Figure 1 1 is a flowchart represent- 
ing the processes in omega.f, specifically, the approximations to the Chapman-Enskog method and the 
Gauss-Laguerre quadrature loop. 

The electron momentum-transfer cross sections for the 30 species are given as a function of 
incident electron energy, in units of electronvolts. 11 ’ 12 Using the Gauss-Laguerre quadrature for a given 
species, which uses abscissa and weight factors, represented by and “m” loops, respectively, in 
figure 11, omega.f creates a table of nine “omegas” for each of the 801 temperatures of interest (from 
1,000 to 9,000 K at 10 K intervals). The “omegas” physically represent the different / and 5 values of the 
average collision cross sections from equation (32). The program creates a separate file for each species, 
containing that specie’s collision cross-section values, which are later used in the calculation of the q mn 
elements, such as in equation (31). For each specie cross-sectional table, linear interpolation is used to 
obtain cross sections for the parameters of the actual potentials (i.e., the actual temperatures of inter- 
est). The overall electrical conductivity calculation procedure is as follows: (1) omega.f calculates the 
“omega” integrals for electron-neutral collisions and creates the omega tables, (2) ceb.f calculates the 
thermodynamic data (e.g., temperature) and species concentrations (including ions), (3) econ.h calcu- 
lates the “omega” integrals for charged-particle collisions and interpolates electron-neutral collision 
cross sections based on temperatures from ceb.f, and (4) econ.h calculates the total electrical conductiv- 
ity, based on equation (25). 

3.2.3 Thermodynamic and Electrical Property Model Implementation 

The purpose of the UMM thermoelectric routines is to provide thermodynamic property data 
(temperature, molecular weight, the ratio of specific heats, laminar viscosity, and total charge) and elec- 
trical transport property data (electrical conductivity data) to the MHD numerical model. The thermo- 
electric codes accomplish this by creating a table that loops in pressure and enthalpy and contains 
the required thermoelectric data, as shown in figure 12. 

As discussed in previous sections, the user specifies in thermhd.f the maximum and minimum 
pressure and enthalpy values, and the number of steps the code should take for the pressure and enthalpy 
loops. (Note that the code calculates each loop’s step size by taking the difference between the maxi- 
mum and minimum loop values and dividing that by the number of steps that should be taken. By mak- 
ing the step size a calculated value, the user is guaranteed that the thermoelectric codes will use their 
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Species Loop 
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Define Incident Electron Energy 
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Transfer Cross Sections 
for Each Species 

I 

Iterate Cross Sections Based 



Figure 1 1 . Flowchart of omega.f routine for execution of the Gauss-Laguerre quadrature procedure. 


exact minimum and maximum values, and that the step sizes will be constant.) As equation (13) illus- 
trates, for the first pressure step, the thermoelectric codes calculate the required properties at every step 
in the enthalpy loop. Then, the code advances to the second pressure step, calculates all the required 
values at every step in the enthalpy loop, again, and continues to advance in the same manner until both 
the enthalpy and pressure loops are completed. Currently, the maximum number of steps for pressure 
and enthalpy are 60 and 500, respectively; however, both values can be increased if necessary. 
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Figure 12. Structure of the thermoelectric property table created by the thermodynamic 
and electrical transport property submodels in the UMM code. 


As previously stated, thermhd.f is the “managing” program of the thermoelectric routines. 

It initiates the thermodynamic and/or electrical transport property calculations, runs the main pressure 
and enthalpy loops, and controls most of the important I/O, including writing the thermoelectric data 
output file (illustrated in fig. 12) as required by the MHD block. Figure 13 is a flowchart representing 
the processes in thermhd.f, which includes the calculations made by econ.h and ceb.f. 

3.3 Universal MGMHD Model Modifications and Improvements 

This section details the modifications and additions made to the baseline MGMHD code. 

The resulting modified code, called mhd.f (with input files mhd.in, mhncom.h, and angle.in) defines 
the MHD block of the UMM code. 

3.3.1 Accelerator Mode 

For a given MHD channel, the primary difference between generator and accelerator modes is 
the direction of the J y current. The reversal of the J y current from generator mode to accelerator mode 
is typically accomplished with an applied current. That is, an external power supply is used to over- 
power the naturally occurring J y current, which stems from the uxB Lorentz force. This applied J y 
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Figure 13. Flowchart for the thermhd.f program — includes the thermodynamic calculations 
made by ceb.f and the electrical calculations made by econ.h. 


current results in a J y B Lorentz force that accelerates the flow in the positive u direction. Furthermore, 
in an MHD accelerator, since the user applies/specifies the J y current instead of one being induced by 
the interaction of the moving particles and the magnetic field, it can alter the solution procedure for the 
problem. 
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The following sections explain how the MGMHD code was modified to achieve MHD accel- 
erator operation in the respective configurations (Faraday, diagonal, and Hall). In all cases, even where 
potentials and currents are calculated in different ways, these modifications were additions to the 
MGMHD code and none of the original MGMHD functionality was lost. It is also important to note that 
when running the UMM code in any MHD accelerator configuration, the user should always double- 
check the output to ensure the proper orientation of the J y current. In accelerator mode, the J y current, 
under certain circumstances, sometimes reverses direction, essentially changing the accelerator to gen- 
erator mode. 


The primary goal of this research is to numerically model a diagonal MHD accelerator. There- 
fore, the majority of the modifications to the MGHMD code deal with facilitating the simulation of a 
diagonal accelerator. This section will focus on the modifications required to convert from a diagonal 
generator model to a diagonal accelerator model. For details of the other modifications and additions, 
including varying of the diagonal angle, axial current neutralized mode, power-takeoff, electrical effi- 
ciency, and three-dimensional total pressure, see the detailed descriptions in Turner. 10 

It is important to note that the governing equations and assumptions outlined in the “Baseline 
Modeling Capabilities” continue to hold. As explained in figures 9 and 10, the physical difference 
between a diagonal accelerator and generator is the diagonal link between the electrodes. Therefore, 
for an MHD accelerator, a new set of oblique coordinates must be defined, and for the purposes of this 
derivation, a new diagonal angle, a d , will be used. Using these new oblique coordinates and a d , it is 
possible to define the electric potential as 


\j/=-(x+ycota d )E x +(p(y,z) , (34) 

where the relationship a d and 6 d is as follows: 

a d =m~0 d . (35) 

Note that equation (12) is recovered when equation (35) is substituted into equation (34). Sub- 
stituting equation (34) into equations (8)— (10) results in the following equation for the two-dimensional 
potential function: 


dy 


do 
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+ dz 






d(p 

dz 


= E xfy[(P + ^ a d)°n]-^[ a n B ( u ~H] ’ 


dy 1 


(36) 


which differs from equation (13) only in the angle and the sign preceding the cotangent. Therefore, 
the same linear decomposition used on equation (13) can be applied to equation (36). 


In the accelerator case as in the generator case, both <p x and (p 2 are independent of E x , which is 
determined by the external connection and loading of the accelerator. Applying Ohm’s Law, the current 
densities are as follows: 


J =E (7 
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(37) 
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where, again, the only differences from the generator equations are the angle and sign before 
the cotangent. The cross-sectional, spatial-averaged current densities are represented as follows: 


(39) 



For the boundary conditions of a diagonal connection, the axial electric field is determined 
by the following equation: 


J 


Id 



(P ~ COta d)(<P2y + uB )) ~ B { (7 n f 1 + P COta d) v ) 
(l + cot 2 a d )(<7 W (/? - cota^ )*P\y) 


(42) 


The purpose of this brief derivation was to show the differences and similarities in the generator 
and accelerator cases. A quick comparison of equations (18) and (37), equations (19) and (38), equations 
(22) and (40), equations (23) and (41), and equations (24) and (42), respectively, shows that the only 
difference between the generator and accelerator cases is the sign preceding the cotangent of the angle. 
Here, it is important to recall a basic property of the cotangent. That is, the cotangent of an angle is equal 
to the negative value of the cotangent of the supplement of that angle. In other words, if one uses an 
obtuse diagonalization angle (if 6 d is >180°), the result would be the same as if one used the angle a d , 
described in this section. Therefore, in order to modify the MGMHD code to run in accelerator mode, 
one needs only to define the diagonalization angle in terms of the obtuse angle, 6 d . Thus, in the MHD 
block of the UMM numerical model (i.e., mhd.f), the baseline MGMHD code was altered to accept 
obtuse diagonalization angles, which the user inputs in mhd.in. 


3.3.2 Integration of Thermodynamic and Electrical Property Models 

Figure 13 shows the structure of the UMM thermoelectrical data table, while figure 14 illustrates 
how these thermoelectrical data are used by the MHD code in three-dimensional calculations. The UMM 
MHD block solves for pressure at each axial location using the value of static pressure at the previous 
axial station and the pressure gradient from the momentum equation, and solves for enthalpy in each 
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Figure 14. Illustration of UMM MHD block using the thermoelectrical data 
in three-dimensional calculations. 


cell at each axial position using the energy equation. (See sec. 3.1.1, Fluid Dynamic Model). Similarly, 
using the momentum equation and mass continuity equation, the UMM MHD block solves for the value 
of axial velocity in each cell at each axial location. For each cell at each axial location, the code inter- 
polates for the values of temperature (T), molecular weight {MW), ratio of specific heats (y), laminar 
viscosity (/q), electrical conductivity (cr), and total charge (eN e ) using the calculated values of pressure 
and enthalpy and the values of pressure and enthalpy from the thermoelectrical data tabulation. With this 
information and the aforementioned axial velocity measurements, virtually any required magnetohydro- 
dynamic parameter (e.g., gas constant, sonic speed, Mach number, total pressure, and Hall parameter) 
may be calculated. Within the mhd.f program, tabular data look-up, two-dimensional interpolation, 
and pointwise calculation of thermoelectrical properties in the cross-sectional plane occur in different 
subroutines. 
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4. MAGNETOHYDRODYNAMIC AUGMENTED PROPULSION EXPERIMENT 
REGIME AND ANALYSIS PARAMETERS 


Because the point of focus of this research is directed at three-dimensional simulation of the 
MAPX accelerator, all other external factors are considered frozen for this analysis. In other words, all 
upstream conditions are considered to be constant as defined at the accelerator channel inlet. There- 
fore, all user-defined input parameters to the thermoelectric codes will remain constant in the follow- 
ing MAPX analyses. The pressure range, enthalpy range, working fluids, weight percentages, and inlet 
temperature are all predetermined by the arc-heater stagnation conditions and primary nozzle expan- 
sion characteristics. Therefore, the same thermoelectric data table was used for all MAPX performance 
analyses. This is also helpful in determining the causes of certain phenomena observed in the accelera- 
tor channel calculations. After a synopsis of the input and flow parameters, the results of an exploratory 
analysis over the trade-space were presented followed by detailed three-dimensional calculations for 
the selected design point. 


4.1 Input and Flow Parameters 

The thermodynamic flow parameters are presented in table 3. These values are input for the 
thermodynamic models (the modified CEA code, called ceb.f, and the managing thermodynamic file, 
thermhd.f). All of these values are held constant throughout the analysis described in this TM. 


Table 3. MAPX defined thermodynamic parameters. 


Working fluid, N 2 
Working fluid temperature 
Seed, Na 
Seed, K 

Seed, NaK temperature 
Seeding percentage 
Pressure range 
Pressure steps 
Enthalpy range 
Enthalpy steps 


100 wt.% 

2,700 K 
22 wt.% 

78 wt.% 

300 K 
1 % 

0.001-5 bar 
60 

7.0x10 5 - 3.5x1 0 7 J/kg 
500 


Table 4 shows the entrance input parameters as they are in the MHD input file. These parameters 
define the flow that enters the MAPX accelerator, and define some of the physical parameters and con- 
nections of the accelerator. For all analyses presented in this TM, these values were held constant. Note 
that “PTO steps” defines the number of staircase steps in the power-takeoff current distribution at the 
entrance of the MAPX channel. 


29 


Table 4. MAPX defined input parameters to UMM. 


Flow rate 
Flow temperature 
Wall temperature 
Static pressure 
Velocity in y-direction 
Velocity in z-direction 
Wall roughness height 
Axial pressure gradient 
Electrode pitch 

PTO length (% of total length) 
PTO steps 


0.130 kg/s 
2,700 K 
1,000 K 

3.24240x1 0 5 Pa 


0.015 m 
0.083% 
5 


Figure 15 shows the physical characteristics of the MAPX channel. These are the same values 
presented in table 1; however, figure 15 depicts the difference in the length between the actual powered 
length and the total channel length. Note that the powered length is 9 cm shorter than the total channel 
length, and the UMM code only models the powered region of the channel. It is also important to note 
that the MAPX channel begins and ends with an insulator. The magnetic field used in these analyses 
is presented in figure 7. 


Electrodes 



Not to scale. 
Dimensions in cm. 


Figure 15. Physical characteristics of the MAPX accelerator. 


4.2 Thermodynamic/Electrophysical Regime 

Figure 16 shows some of the major parameters over the regime of interest to the MAPX accel- 
erator. Please note that in order to see details in each plot, the axes are not oriented in the same direction. 
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Figure 16. Major thermodynamic/thermoelectric parameters over regime of interest: (a) Temperature 
as a function of pressure and enthalpy, (b) electrical conductivity as a function of tempera- 
ture and pressure, (c) electron number density as a function of temperature and pressure, 
and (d) electron mobility as a function of temperature and pressure. 


From figure 16, temperature increases with increasing enthalpy, but seems to remain nearly 
constant for a given enthalpy with increasing pressure. The same is true for electrical conductivity — as 
temperature increases, conductivity also increases (after temperature reaches « 2,500 K). However, as 
pressure increases, this seems to cause the conductivity to decrease slightly, for a given temperature. 
This is logical — increasing pressure would inhibit electron movement and therefore conductivity. In 
part (c), electron number density increases with increasing temperature (after 3,000 K) and pressure. 
This is also logical — increasing temperature releases more electrons, and increasing pressure adds more 
particles per unit area. Electron mobility also seems to increase with temperature, but decrease with 
increasing pressure. 
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5. PARAMETRIC TRADE STUDY AND DESIGN OPTIMIZATION ANALYSES 


In the parametric analysis of the trade space, applied currents from zero through 300 A (25- A 
increments) were tested against different diagonal angles ranging from 0° to 60°, as measured from the 
positive y-axis. In addition to this diagonal angle trade study, the UMM code was also run in “axial cur- 
rent neutralized mode,” which, in principle, should provide a reliable indication of the optimum angle 
for superior performance given a specific channel configuration. Since only the J y B Lorentz force pushes 
the flow downstream, axial current provides no positive acceleration effect and therefore reduces overall 
performance. Hence, the axial current neutralized mode, which uses the cross-sectional averaged current 
values, yields an estimate of the optimum diagonalization angle for best performance. 

Figure 17 shows the axial variation of diagonalization angle, 6 d , for axial current neutralized 
case with 100-A applied current. According to this figure, the optimum diagonalization angle (d d , as 
measured from the positive x-axis) is between 140° and 150° over most of the channel. However, due 
to physical constraints, a 135° diagonalization angle is the maximum that can be practically attained 
in the accelerator, and is therefore chosen for the remainder of the parametric analyses in this TM. 

It is important to how the diagonalization angle is defined. In this case, 6 d is measured from 
the positive x-axis and is always obtuse in accelerator mode. However, the diagonalization configuration 
may also be adequately defined with reference to the electric field angle, 6, which is measured from the 
positive y-axis. This angle is always an acute value (for generators and accelerators), and one must sim- 
ply know whether a generator or accelerator is being discussed. Therefore, the 135 diagonalization angle 
(8 d ) is equivalent to a 45 electric field angle (0). 

Figure 18 shows the plots of cross-sectional averaged (a) velocity and (b) temperature as a func- 
tion of downstream axial location in the MAPX accelerator. The different lines and numbers on each 
plot represent different levels of applied current (zero to 300 A, where the “hydro” case has no applied 
current and no magnetic field). From figure 18(a), it is observed that the highest levels of applied cur- 
rent do not result in the highest cross-sectional averaged axial velocities. This is due, in large part, to 
the excessive joule heating under high current loading conditions, as indicated in figure 18(b). Specifi- 
cally, the highest cross-sectional averaged velocity occurs in the 125- A case. Higher current levels show 
reduced axial velocities. Part (b) shows that both the 125- and 100-A applied current cases have very 
stable cross-sectional temperatures. After some initial heating, they remain constant through the remain- 
der of the channel. The highest levels of applied current show excessive and runaway heating and 
detrimental performance effects. Also note that in some cases (e.g., 300 A), the UMM code terminated 
before reaching the end of the channel due to the occurrence of flow reversal, which violates the PNS 
approximation. 

Figure 19 shows the centerline total pressure as a function of the downstream axial distance at 
various levels of applied current. The flow enters the channel with a total pressure of ~8.4 atm. Depend- 
ing on the level of applied current, the flow exits the accelerator channel with a net gain or loss in total 
pressure. 
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Figure 17. Axial variation of diagonalization angle, Q d , for axial current neutralized case 
with 100-A applied current. Definition of 6 and 0 d illustrated above. 


Observed was that the highest levels of applied current (175 to 300 A) experience the severest 
losses in total pressure. This is consistent with figure 1 8(a) in which the highest levels of applied current 
cause axial velocity deceleration in the second half of the channel. Furthermore, all levels of applied cur- 
rent experience some form of total pressure loss in the channel. In the 100-A case, this total pressure loss 
occurs =75 cm downstream, suggesting that frictional forces begin to overpower the flow at this loca- 
tion. As seen in figure 1 8(a), the slope of the axial velocity curve also begins to flatten out at this loca- 
tion, suggesting that the decrease in axial acceleration is related to the decrease in total pressure. 

The cross-sectional averaged J y current density is shown in figure 20 at various levels of applied 
current. Note that the stair-step distribution at the entrance to the channel is associated with the power- 
takeoff region. Of particular relevance, however, is the fact that J y tends to decrease down the channel. 
This is a critical effect since it is J y that controls the magnitude of the JxB Lorentz force and thereby 
determines acceleration effectiveness. The resulting J current arises from the net electric field, which is 
the difference between the induced electric field (uxB) and the externally applied field. Therefore, as 
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Figure 18. Cross-sectional averaged (a) velocity and (b) temperature as a function of axial 
distance in the MAPX accelerator. Levels of applied current are indicated. 
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Figure 19. Centerline total pressure as a function of axial position. Levels of applied current 
are indicated. 


flow velocity increases (as indicated in fig. 1 8), the induced electric field increases, yielding an overall 
net reduction in the transverse electric field. The net decrease in the transverse electric field therefore 
results in a net decrease in the transverse current, although total current is conserved. It is also important 
to realize that J y is a measure of current density, and as the cross-sectional area of the channel increases, 
the current density will naturally decrease, while the velocity will naturally increase due to supersonic 
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Figure 20. Cross-sectional averaged transverse current density, J y , as a function of axial position. 
Levels of applied current are indicated. 


expansion. Based on the results presented in figure 20, the 100- and 120-A applied current cases main- 
tain the highest levels of J y current density throughout most of the MAPX channel, and are therefore 
near optimum from a performance perspective. 

Figure 21 shows the axial current density, J x , as a function of axial distance — part (a) showing 
the cross-sectional average J x values and part (b) giving the centerline values for J x . The differences 
between the cross-sectional averaged and centerline axial current densities are significant. The cross- 
sectional averaged current densities are all <0.05 A/cm 2 , while the centerline values for J x are as high 
as 70 A/cm 2 . The large centerline values coupled with the very small cross-sectional averaged values 
indicate that within each cross-sectional plane, there exist upstream and downstream axial currents of 
approximately the same magnitude. 

Recall that the axial current neutralized mode uses the cross-sectional averaged values of J x to 
calculate the appropriate angle. Since the diagonal angle is applied across the entire cross section, it is 
appropriate for the algorithm to use the cross-sectional values of J x . The nonzero values for the cross- 
sectional averaged J x , as indicated in figure 21, are a result of the physical limitations of the diagonal 
angle configuration. Because the recommended axial current neutralized angle of 55° is not physically 
practical, a 45° angle was selected for evaluation, hence, the observance of small but finite cross- 
sectional averaged J x values. The axial current neutralized mode does, however, allow for localized 
nonzero values of J x , and the local values can be quite large as shown in figure 21. 

Power levels are also of great importance when designing and optimizing an MHD accelerator. 
Figure 22(a) shows the total cumulative power that the MHD accelerator delivers to flow at various 
levels of applied current. Note that the pure hydrodynamic and unpowered cases show no power deliv- 
ered by the MHD channel. Also note that, for the cases that did not terminate prematurely, there is a 
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Figure 21. Axial current density, J x , as a function of axial position: (a) Cross-sectional averaged 
J x and (b) centerline J x . Levels of applied current are indicated. 


diminishing rate of return for increasing levels of applied current. For example, between 75 and 
100 A, the difference in total cumulative power is 68 kW; whereas, between 100 and 125 A, the differ- 
ence in total cumulative power is =45 kW. As can easily be seen on the plot, the increase in total power 
diminishes for each additional 25-A increment for those cases for which the calculation did not terminate 
before the end of the channel. Since this graph shows a rate of diminishing total power with increased 
applied current, one would assume that there is an optimum level of applied current, above which losses 
would dominate. 

The total cumulative power delivered to the flow, as shown in figure 22(a), either serves to accel- 
erate the flow or heat the flow. Figure 22(b) shows the amount of total power delivered to the flow that 
was used to accelerate, or actually “push” the flow, as accumulated over the length of the MAPX accel- 
erator. This graph shows the same diminishing returns on push power as were seen for total power; how- 
ever, in this case, the 150- and 175-A cases show that push power drops off in the downstream region. 
Note, however, that between 75 and 100 A, the difference in push power is =36 kW, while the difference 
in push power between 100 and 125 A is only =15 kW. So, although the 125-A level of applied current 
does represent the highest level of accumulated push power, the difference from the previous level of 
applied current is less. It is also interesting to note that the unpowered case shows a small, but negative, 
push power. This is consistent with having no applied current, but allowing the magnetic field to induce 
a small, internal negative J y current, even though this is an open-circuit case. 

Aside from acceleration, the total applied power can serve to heat the flow through joule heat- 
ing or joule dissipation. Figure 22(c) shows the amount of the total power that was used to heat the flow 
providing no propulsive benefits whatsoever. The graph shows a similar “diminishing returns” phenom- 
enon as seen in parts (a) and (b); however, the 150- and 175-A cases show enormous, increasing jumps 
in dissipated heat losses. Furthermore, the difference between the 75- and 100-A cases and between the 
100- and 125-A cases is almost negligible. Between 75 and 100 A, the difference is =32 kW, while the 
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Figure 22. Various power levels along the channel: (a) Accumulated total power, 

(b) push power, and (c) joule dissipated power. 
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difference between the 100- and 125-A cases is -30 kW. Also note that the unpowered case shows 
a small, positive joule dissipation loss equal and opposite to the push power in (b), thereby yielding 
a zero summation. 

Figure 22 tells an important story. While it is true that the 125-A case has the highest accumu- 
lated push-power level, one should note that only 1 5 kW of the 45 kW added total power (above the 
100-A case) went toward flow acceleration. The other 30 kW, twice the push-power value, went directly 
into joule dissipation losses (i.e., heating). Thus, only one-third of the added total power between 100 
and 125 A was used to accelerate the flow at the 125-A level. Moreover, -52% of the added total power 
between the 75- and 100-A current levels was used for acceleration in the 100-A case. Although the 
125-A case does give higher total and push-power values, more of the additional applied power (above 
the 100-A case) is used to heat the flow rather than to accelerate the flow. 

The differences in the push and applied power levels can best be explained through the global 
electrical efficiency, shown in figure 23. This graph clearly shows that, from a standpoint of power utili- 
zation, the 100-A case is superior to the 125-A case. Furthermore, since the 75- and 125-A cases are so 
near to each other (and both approximately equidistance from the 100-A case), the 100-A case seems to 
represent the highest electrical efficiency (i.e., 46%) for the MAPX accelerator configuration. It should 
also be noted that the initial decrease in global efficiency in figure 23 is due to the low current levels in 
the power-takeoff region. In short, for an MHD accelerator used as an onboard propulsive device, effi- 
ciency is paramount, and the 100-A level of applied current yields the best efficiency. 



Figure 23. Cumulative global efficiency for the MAPX accelerator. Levels of applied current 
are indicated. 
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6. PRETEST DESIGN POINT ANALYSIS 


The optimal MAPX design point, as predicted by the UMM numerical model, occurs near 100-A 
applied current when diagonally connected with 6 = 45° (0 d = 135°). Under these conditions, the cross- 
sectional averaged velocity should almost double, the centerline total pressure should increase by =33%, 
and the cross-sectional averaged temperature should remain approximately constant following a slight 
initial increase. The overall accelerator efficiency is =46%. This section provides a detailed discussion 
of the three-dimensional simulation for this optimized design point. 

Figure 24 shows the three-dimensional axial velocity at various downstream locations. The 
UMM code treats the flow at the entrance of the MAPX channel like a plug flow; therefore, the plot at 
0.2 cm is still plug shaped. As the flow develops downstream, however, one notices that the bulk of the 
flow is attracted to the bottom wall (i.e., the electron accepting electrode surface). Therefore, the bulk of 
the flow seems to be pushed in the direction of the transverse electron flow. 

Figure 25 shows a contour plot of the axial velocity, with a vector overlay representing the sec- 
ondary flow velocity components (i.e., v and w). MHD effects tend to increase the amount of secondary 
flow in a channel, as evident in figure 25. Generally, the secondary flow seems to be driving the axial 
core flow toward the bottom electrode. The center of the secondary flow, which appears to feed off the 
high-velocity regions in the primary axial flows, moves downward in the channel as the flow develops 
downstream. This behavior has been previously observed in MHD channels. 18 The upward-directed sec- 
ondary flow regions seem to be pushing the axial flow into the upper comers of the cross section, while 
simultaneously pushing the flow away from the center of the top electrode. This cross-plane circula- 
tion causes a concave region to form at the top of the axial profile, which, if allowed to develop further, 
might cause a separation in the core flow field. 

Figure 18 shows that, following an initial increase, the cross-sectional averaged temperature 
profile for the 100-A case remained fairly constant. Figure 26, however, tells a different story. Figure 26 
shows that, while the cross-sectional averaged temperature may be constant, there are large fluctuations 
in the temperature profile as the flow moves downstream. Most noticeable is the hot spot that begins to 
develop along the bottom electrode (i.e., the electron accepting electrode surface). This high-temperature 
region also corresponds with the downward skewing of the velocity profile and the strong secondary 
flow effects. 

Figure 19 demonstrated that the centerline total pressure increased by =30% above the inlet con- 
dition. However, figure 27 shows that, while centerline total pressure does increase, there are severe vis- 
cous losses near the channel sidewalls. The flow is forced to enter as a plug, but sidewall shear quickly 
acts to retard the near-wall flow, thereby degrading total pressure in the region. Near the exit of the chan- 
nel, the viscous losses from the sidewalls have actually started to encroach into the central region of the 
channel, limiting total pressure gains to a small region near the centerline. 
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Figure 24. Three-dimensional axial velocity, u, plotted at (a) 0.2, (b) 10, (c) 25, (d) 60, (e)80, 
and (f) 90.5 cm along the MAPX accelerator channel. The last position corresponds 
to the end of the powered region in the channel. The positive J y current density flows 
from the bottom electrode (background) to the top electrode (foreground) in a direction 
out of the page. All velocities are in m/s. 
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Axial Velocity (m/s) Axial Velocity (m/s) 






u: 179 357 536 714 893 1,071 1,250 1,428 1,607 1,786 1,964 2,143 2,321 2,500 2,678 


Figure 25. Vector plots of cross-sectional secondary flow velocities for (a) 0.2 cm, (b) 10 cm, 
(c) 25 cm, (d) 60 cm, (e) 80 cm, and (f) 90.5 cm with the axial velocity, u, shown 
as a background contour plot. Positive J current flows from the bottom electrode 
to the top electrode. The length of the vectors is representative of their relative 
intensities. All velocities are in m/s. 
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Figure 26. Three-dimensional temperature plots at various downstream locations in the MAPX 
accelerator channel: (a) 0.2 cm, (b) 10 cm, (c) 25 cm, (d) 60 cm, (e) 80 cm, and 
(f) 90.5 cm. Positive J y current density flows from the bottom electrode (background) 
to the top electrode (foreground) in accelerator mode. 
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Figure 27. Three-dimensional total pressure plots at various downstream locations: (a) 0.2 cm, 
(b) 10 cm, (c) 25 cm, (d) 60 cm, (e) 80 cm, and (f) 90.5 cm. Positive J y current flow 
is from the background to the foreground. 
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Figure 20 shows that the cross-sectional averaged transverse current density, J y , declines along 
the accelerator channel, but never turns negative. Figure 28 reiterates that point with a contour plot of 
the transverse current density. The most intense current density occurs at =10 cm. Beyond that axial 
location, the current density declines, and the lowest current densities occur near the center of the chan- 
nel. Again, a portion of the decline in current density may be simply attributed to the increase in channel 
cross-sectional area. 

Because J y and J z currents flow in the cross-sectional plane, they can be accurately visualized in 
relation to the secondary fluid flow. Figure 29 shows the total current flow in the cross-sectional plane, 
which is basically the vector sum of the J y and J z current densities. The most important aspect of the 
cross-sectional current is that it does not reverse direction and begin flowing in the negative direction. 

As long as J y is positive, the J y B Lorentz force still provides some push to the flow. If the current 
reverses direction to flow negative, then the Lorentz force would act to decelerate the flow. 

In comparison to figure 29, figure 30 shows the vector plot of cross-sectional current densities 
for the 175 -A applied current level. Here, the high applied current level yields a positive J y somewhat 
beyond the 40-cm location. At the 60-cm position, however, J y begins flowing in the negative direction 
and continues to do so through the 80-cm location. This phenomenon sets up a decelerating force in 
opposition to the bulk flow direction. 

One explanation for the current reversal depicted in figure 30 is the intense heat generation 
shown in figure 3 1 . Starting as early as 45 cm, the lower half of the channel begins to experience exces- 
sive heating. Between 70 and 90 cm, an intense hot spot develops along the bottom electrode, much 
like the 100-A case of figure 26, but to a lesser degree. The intense heating at 175 A may also be seen 
in figure 18, thereby explaining the low efficiency, current reversal, and overall poor performance. 

Figure 21 shows the centerline and cross-sectional averaged axial current density, J x , as a func- 
tion of axial location. The large differences between the centerline and cross-sectional averaged current 
densities in figure 21 indicate strong localized, opposing currents in the x-direction. Figure 32 confirms 
that the MAPX accelerator channel does have very intense opposing local axial currents in the cross- 
sectional plane. At the 80-cm position, for instance, very strong, local positive currents exist along the 
bottom electrode. 
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Figure 28. Contour plots of the J y current density (A/cm 2 ) at various downstream locations: 
(a) 0.2 cm, (b) 10 cm, (c) 25 cm, (d) 60 cm, (e) 80 cm, and (f) 90.5 cm. Positive J y 
current flow goes from the bottom electrode to the top electrode. 
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Figure 29. Vector plots of cross-sectional current densities J y and J z at various downstream 
locations in the MAPX accelerator channel for the 100- A case: (a) 0.2 cm, 

(b) 10 cm, (c) 25 cm, (d) 60 cm, (e) 80 cm, and (f) 90.5 cm. Positive J y current 
flow is from the bottom to top electrodes. 
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Figure 30. Vector plots of the cross-sectional current flow at various downstream locations 
for the 175-A case: (a) 45 cm, (b) 60 cm, (c) 65 cm, (d) 70 cm, (e) 80 cm, and 
(f) 90.5 cm. The positive J y current density is from the bottom to top electrode 
when in accelerator mode. 
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Figure 3 1 . Contour plot of channel temperature (K) at various downstream locations 
for the 175-A case: (a) 25 cm, (b) 45 cm, (c) 60 cm, (d) 70 cm, (e) 80 cm, 
and (f) 90.5 cm. Positive J y current flow occurs from the bottom to top 
electrode when in accelerator mode. 


48 








(a) 0.2 cm 


*% 

(b) 10 cm 


% 

(c) 25 cm 



Figure 32. Three-dimensional axial current, J x , plots at various downstream locations: 
(a) 0.2 cm, (b) 10 cm, (c) 25 cm, (d) 60 cm, (e) 80 cm, and (f) 90.5 cm. 

The positive J y current flows from the bottom electrode (background) to top 
electrode (foreground) in accelerator mode. 
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7. CONCLUSION 


A three-dimensional numerical MHD model was developed for the purpose of analyzing and 
optimizing the MAPX. The results of this analysis were intended to offer a better understanding of 
the critical physical processes within the MHD accelerator channel, to assist in the identification of an 
optimized configuration for initial testing, and to provide detailed pretest performance predictions of that 
optimized configuration. 

The three-dimensional nature of the coupled fluid dynamic and electromagnetic effects in MHD 
flows, particularly enclosed channel flows, requires the use of a proper three-dimensional numerical 
model. Unfortunately, no such model was generally available for accelerator configurations at the begin- 
ning of the MAPX project. Therefore, in order to achieve the capabilities and accuracy required for a 
comprehensive analysis, modifications and additions were made to an existing three-dimensional PNS 
model previously developed by the Argonne National Laboratory for MHD generator channels. The 
resulting model, with integrated thermodynamic and electrical transport property routines established on 
the basis of a rigorous kinetic theory treatment with reliable collisional cross-section data, was then used 
to simulate the partially ionized flow inside the MAPX accelerator channel. 

The results of an extensive three-dimensional analysis and trade study of the MAPX accelerator 
revealed critical effects affecting performance and provided assistance in the definition of an optimized 
configuration for initial testing. These results clearly indicate that entropy-producing mechanisms domi- 
nate the flow due to the high surface-to-volume ratio of the MAPX channel. More specifically, accel- 
erator performance was found to be extremely sensitive to turbulent shear losses, which are difficult 
to simulate with high precision. Furthermore, large local axial currents cause asymmetries in the flow, 
resulting in localized runaway heating. The results from this extensive performance optimization study 
indicated that the MHD accelerator should utilize a 45° diagonalization angle with the applied current 
evenly distributed over the first five inlet electrode pairs. When powered at 100 A, this configuration is 
expected to yield a 50% global efficiency with an 80% increase in axial velocity and a 50% increase in 
centerline total pressure. 
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